###R code US electricity###

#front matter
rm(list=ls())
#library(TSA)
library(foreign)
library(lattice)
library(plm)
library(stargazer)
library(xtable)
#library(XLConnect)
library(reshape)

#set working directory
#setwd("E:/stateCO2")
#setwd("/Volumes/MONOGAN/stateCO2")

#Cleaning the data
#gdp.wide<-read.csv("Archive1/gdpWide.csv")
#annual.vars<-names(gdp.wide)[-1]
#gdp.long<-melt(gdp.wide,id.vars="GeoName",measure.vars=annual.vars)
#write.csv(gdp.long, "Archive1/gdpLong.csv",row.names=F)
#gdp.long<-read.csv("Archive1/gdpLong.csv")
#year<-c(1990:2014)
#multiplier<-c(.561,.585,.602,.620,.636,.654,.674,.689,.700,.715,.739,.760,.772,.790,.811,.838,.865,.890,.924,.921,.936,.966,.986,1.000,1.016)
#inflation<-cbind(year,multiplier)
#gdp.long<-merge(x=gdp.long,y=inflation,by="year")
#gdp.long$realGDP<-gdp.long$nominalGDP/gdp.long$multiplier
#library(datasets); data(state)
#gdp.long$state<-state.abb[match(gdp.long$GeoName,state.name)]
#gdp.long[1,]
##gdp.long[order(gdp.long$realGDP),]
#POLICYII<-read.csv("UsElectric9.csv")
#POLICYII<-merge(x=gdp.long,y=POLICYII,by=c("state","year"))
#POLICYII$gdp<-log(POLICYII$realGDP)
#POLICYII$ee<-POLICYII$sale/POLICYII$gdp
#write.csv(POLICYII,"UsElectric10.csv",row.names=F)

#load data
POLICYII <-read.csv("UsElectric.csv") 

#create two more variables
POLICYII$rpsprop<-POLICYII$rpsdata/exp(POLICYII$gdp)
POLICYII$co2gdp<-exp(POLICYII$pop)*POLICYII$co2pop/exp(POLICYII$gdp)

#table of descriptives
used.vars<-subset(POLICYII,select=c(co2gdp,renew2,ee,dere,ipp,rpsprop,isp,gdp,price,gasPrice))
desc<-cbind(
apply(used.vars,2,mean),
apply(used.vars,2,sd),
apply(used.vars,2,min),
apply(used.vars,2,max)
)
rownames(desc)<-c("CO2 per GDP","Proportion renewable","Electricity per GDP","Retail Deregulation","Proportion independent","Renewable portfolio standard","RTO or ISO present","Logged state GDP","Electricity price","Gas price")
colnames(desc)<-c("Mean","Std. Dev.","Minimum","Maximum")
print(xtable(desc,caption="Descriptive Statistics, 50 states from 1990-2014",digits=4),"HTML",caption.placement="top")

#quick look at states' policies
table(POLICYII$dere,POLICYII$state)
table(as.numeric(POLICYII$rpsprop>0),POLICYII$state)

###MODELS###
mod1=plm(price~ rpsprop+dere+ipp+isp+gasPrice,data=POLICYII,index=c("state","year"),model="random");summary(mod1)
mod2=plm(qnorm(renew2+.00001)~rpsprop+dere+ipp+isp+gdp+gasPrice,data=POLICYII,index=c("state","year"),model="random");summary(mod2)
mod3=plm(ee~ dere+ipp+isp+price+gasPrice,data=POLICYII,index=c("state","year"),model="random");summary(mod3)
mod4=plm(log(co2gdp)~ rpsprop+dere+ipp+isp+gdp+renew2+ee+gasPrice,data=POLICYII,index=c("state","year"),model="random");summary(mod4)
mod4alt=plm(log(co2pop)~ rpsprop+dere+ipp+isp+gdp+renew2+ee+gasPrice,data=POLICYII,index=c("state","year"),model="random");summary(mod4)

#output table
stargazer(mod1,mod2,mod3,mod4,title="Recursive Model of Renewables, Electricity Consumption, and CO2 Emissions in the 50 States, 1990-2014",star.cutoffs=.1,column.labels=c("Price of electricity","Probit of proportion renewable","Electricity per GDP","Logged CO$^2$ per GDP"),covariate.labels=c("Renewable portfolio standard","Deregulation","Proportion independent","RTO or ISO present","Logged state GDP","Electricity price","Proportion renewable","Electricity per GDP","Gas price","Constant"),omit.stat=c("f","adj.rsq","n"),notes="N=50 and T=25 for 1,250 state-years. Estimated in R 3.4.1.",notes.align='l',type='html',no.space=T)#,type="latex"

#effect of isp (through proportion renewable, electricity per gdp [itself direct and through price], and direct effect)
no.isp.price<-mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*0+mod1$coef[4]*mean(POLICYII$ipp)+mod1$coef[5]*0+mod1$coef[6]*mean(POLICYII$gasPrice);no.isp.price
no.isp.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));no.isp.renew
no.isp.eGDP<-mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*mean(POLICYII$ipp)+mod3$coef[4]*0+mod3$coef[5]*no.isp.price+mod3$coef[6]*mean(POLICYII$gasPrice);no.isp.eGDP
no.isp<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]*no.isp.renew+mod4$coef[8]*no.isp.eGDP+mod4$coef[9]*mean(POLICYII$gasPrice))

yes.isp.price<- mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*0+mod1$coef[4]*mean(POLICYII$ipp)+mod1$coef[5]*1+mod1$coef[6]*mean(POLICYII$gasPrice);yes.isp.price
yes.isp.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*1+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));yes.isp.renew
yes.isp.eGDP<- mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*mean(POLICYII$ipp)+mod3$coef[4]*0+mod3$coef[5]*yes.isp.price+mod3$coef[6]*mean(POLICYII$gasPrice);yes.isp.eGDP
yes.isp<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*1+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]*yes.isp.renew+mod4$coef[8]*yes.isp.eGDP+mod4$coef[9]*mean(POLICYII$gasPrice))

no.isp;yes.isp
100*(no.isp-yes.isp)/no.isp

#effect of retail deregulation
no.retail.price<- mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*0+mod1$coef[4]*mean(POLICYII$ipp)+mod1$coef[5]*0+mod1$coef[6]*mean(POLICYII$gasPrice);no.retail.price
no.retail.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));no.retail.renew
no.retail.eGDP<- mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*mean(POLICYII$ipp)+mod3$coef[4]*0+mod3$coef[5]*no.retail.price+mod3$coef[6]*mean(POLICYII$gasPrice);no.retail.eGDP
no.retail<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]*no.retail.renew+mod4$coef[8]*no.retail.eGDP+mod4$coef[9]*mean(POLICYII$gasPrice))

yes.retail.price<- mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*1+mod1$coef[4]*mean(POLICYII$ipp)+mod1$coef[5]*0+mod1$coef[6]*mean(POLICYII$gasPrice);yes.retail.price
yes.retail.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*1+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));yes.retail.renew
yes.retail.eGDP<- mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*mean(POLICYII$ipp)+mod3$coef[4]*0+mod3$coef[5]*yes.retail.price+mod3$coef[6]*mean(POLICYII$gasPrice);yes.retail.eGDP
yes.retail<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*1+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]*yes.retail.renew+mod4$coef[8]*yes.retail.eGDP+mod4$coef[9]*mean(POLICYII$gasPrice))

no.retail;yes.retail
100*(no.retail-yes.retail)/no.retail

#effect of proportion independent
mean.ipp.price<- mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*0+mod1$coef[4]*mean(POLICYII$ipp)+mod1$coef[5]*0+mod1$coef[6]*mean(POLICYII$gasPrice);mean.ipp.price
mean.ipp.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));mean.ipp.renew
mean.ipp.eGDP<- mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*mean(POLICYII$ipp)+mod3$coef[4]*0+mod3$coef[5]*mean.ipp.price+mod3$coef[6]*mean(POLICYII$gasPrice);mean.ipp.eGDP
mean.ipp<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]* mean.ipp.renew +mod4$coef[8]* mean.ipp.eGDP +mod4$coef[9]*mean(POLICYII$gasPrice))

sd.ipp.price<- mod1$coef[1]+mod1$coef[2]*mean(POLICYII$rpsprop)+mod1$coef[3]*0+mod1$coef[4]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+mod1$coef[5]*0+mod1$coef[6]*mean(POLICYII$gasPrice);sd.ipp.price
sd.ipp.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));sd.ipp.renew
sd.ipp.eGDP<- mod3$coef[1]+mod3$coef[2]*0+mod3$coef[3]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+mod3$coef[4]*0+mod3$coef[5]*sd.ipp.price+mod3$coef[6]*mean(POLICYII$gasPrice);sd.ipp.eGDP
sd.ipp<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]* sd.ipp.renew +mod4$coef[8]* sd.ipp.eGDP +mod4$coef[9]*mean(POLICYII$gasPrice))

mean.ipp;sd.ipp
100*(mean.ipp-sd.ipp)/mean.ipp

#effect of rps
mean.rps.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*mean(POLICYII$rpsprop)+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));mean.rps.renew
mean.rps<-exp(mod4$coef[1]+mod4$coef[2]*mean(POLICYII$rpsprop)+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]* mean.rps.renew +mod4$coef[8]*mean(POLICYII$ee) +mod4$coef[9]*mean(POLICYII$gasPrice))

sd.rps.renew<-pnorm(mod2$coef[1]+mod2$coef[2]*(mean(POLICYII$rpsprop)+sd(POLICYII$rpsprop))+mod2$coef[3]*0+mod2$coef[4]*mean(POLICYII$ipp)+mod2$coef[5]*0+mod2$coef[6]*mean(POLICYII$gdp)+mod2$coef[7]*mean(POLICYII$gasPrice));sd.rps.renew
sd.rps<-exp(mod4$coef[1]+mod4$coef[2]*(mean(POLICYII$rpsprop)+sd(POLICYII$rpsprop))+mod4$coef[3]*0+mod4$coef[4]*mean(POLICYII$ipp)+mod4$coef[5]*0+mod4$coef[6]*mean(POLICYII$gdp)+mod4$coef[7]*sd.rps.renew +mod4$coef[8]*mean(POLICYII$ee) +mod4$coef[9]*mean(POLICYII$gasPrice))

mean.rps;sd.rps
100*(mean.rps-sd.rps)/sd.rps

###FIXED EFFECTS###
fe1=plm(price~ rpsprop+dere+ipp+isp+gasPrice,data=POLICYII,index=c("state","year"),model="within");summary(fe1)
fe2=plm(qnorm(renew2+.00001)~rpsprop+dere+ipp+isp+gdp+gasPrice,data=POLICYII,index=c("state","year"),model="within");summary(fe2)
fe3=plm(ee~ dere+ipp+isp+price+gasPrice,data=POLICYII,index=c("state","year"),model="within");summary(fe3)
fe4=plm(log(co2gdp)~ rpsprop+dere+ipp+isp+gdp+renew2+ee+gasPrice,data=POLICYII,index=c("state","year"),model="within");summary(fe4)
fe4alt=plm(log(co2pop)~ rpsprop+dere+ipp+isp+gdp+renew2+ee+gasPrice,data=POLICYII,index=c("state","year"),model="within");summary(fe4alt)

#output table
stargazer(fe1,fe2,fe3,fe4,title="Fixed Effects Recursive Model of Renewables, Electricity Consumption, and CO2 Emissions in the 50 States, 1990-2014",star.cutoffs=.1,column.labels=c("Price of electricity","Probit of proportion renewable","Electricity per GDP","Logged CO$^2$ per GDP"),covariate.labels=c("Renewable portfolio standard","Deregulation","Proportion independent","RTO or ISO present","Logged state GDP","Electricity price","Proportion renewable","Electricity per GDP","Gas price"),omit.stat=c("f","adj.rsq","n"),notes="N=50 and T=25 for 1,250 state-years. Estimated in R 3.4.1.",notes.align='l',type='html',no.space=T)#,type="latex"

#effect of isp (through proportion renewable, electricity per gdp [itself direct and through price], and direct effect)
no.isp.price<-fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*0+fe1$coef[3]*mean(POLICYII$ipp)+fe1$coef[4]*0+fe1$coef[5]*mean(POLICYII$gasPrice);no.isp.price
no.isp.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));no.isp.renew
no.isp.eGDP<-fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*mean(POLICYII$ipp)+fe3$coef[3]*0+fe3$coef[4]*no.isp.price+fe3$coef[5]*mean(POLICYII$gasPrice);no.isp.eGDP
no.isp<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]*no.isp.renew+fe4$coef[7]*no.isp.eGDP+fe4$coef[8]*mean(POLICYII$gasPrice))

yes.isp.price<- fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*0+fe1$coef[3]*mean(POLICYII$ipp)+fe1$coef[4]*1+fe1$coef[5]*mean(POLICYII$gasPrice);yes.isp.price
yes.isp.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*1+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));yes.isp.renew
yes.isp.eGDP<- fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*mean(POLICYII$ipp)+fe3$coef[3]*0+fe3$coef[4]*yes.isp.price+fe3$coef[5]*mean(POLICYII$gasPrice);yes.isp.eGDP
yes.isp<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*1+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]*yes.isp.renew+fe4$coef[7]*yes.isp.eGDP+fe4$coef[8]*mean(POLICYII$gasPrice))

no.isp;yes.isp
100*(no.isp-yes.isp)/no.isp

#effect of retail deregulation
no.retail.price<- fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*0+fe1$coef[3]*mean(POLICYII$ipp)+fe1$coef[4]*0+fe1$coef[5]*mean(POLICYII$gasPrice);no.retail.price
no.retail.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));no.retail.renew
no.retail.eGDP<- fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*mean(POLICYII$ipp)+fe3$coef[3]*0+fe3$coef[4]*no.retail.price+fe3$coef[5]*mean(POLICYII$gasPrice);no.retail.eGDP
no.retail<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]*no.retail.renew+fe4$coef[7]*no.retail.eGDP+fe4$coef[8]*mean(POLICYII$gasPrice))

yes.retail.price<- fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*1+fe1$coef[3]*mean(POLICYII$ipp)+fe1$coef[4]*0+fe1$coef[5]*mean(POLICYII$gasPrice);yes.retail.price
yes.retail.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*1+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));yes.retail.renew
yes.retail.eGDP<- fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*mean(POLICYII$ipp)+fe3$coef[3]*0+fe3$coef[4]*yes.retail.price+fe3$coef[5]*mean(POLICYII$gasPrice);yes.retail.eGDP
yes.retail<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*1+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]*yes.retail.renew+fe4$coef[7]*yes.retail.eGDP+fe4$coef[8]*mean(POLICYII$gasPrice))

no.retail;yes.retail
100*(no.retail-yes.retail)/no.retail

#effect of proportion independent
mean.ipp.price<- fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*0+fe1$coef[3]*mean(POLICYII$ipp)+fe1$coef[4]*0+fe1$coef[5]*mean(POLICYII$gasPrice);mean.ipp.price
mean.ipp.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));mean.ipp.renew
mean.ipp.eGDP<- fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*mean(POLICYII$ipp)+fe3$coef[3]*0+fe3$coef[4]*mean.ipp.price+fe3$coef[5]*mean(POLICYII$gasPrice);mean.ipp.eGDP
mean.ipp<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]* mean.ipp.renew +fe4$coef[7]* mean.ipp.eGDP +fe4$coef[8]*mean(POLICYII$gasPrice))

sd.ipp.price<- fixef(fe1)["CA"]+fixef(fe1)["CA"]*mean(POLICYII$rpsprop)+fe1$coef[2]*0+fe1$coef[3]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+fe1$coef[4]*0+fe1$coef[5]*mean(POLICYII$gasPrice);sd.ipp.price
sd.ipp.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));sd.ipp.renew
sd.ipp.eGDP<- fixef(fe3)["CA"]+fixef(fe3)["CA"]*0+fe3$coef[2]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+fe3$coef[3]*0+fe3$coef[4]*sd.ipp.price+fe3$coef[5]*mean(POLICYII$gasPrice);sd.ipp.eGDP
sd.ipp<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*(mean(POLICYII$ipp)+sd(POLICYII$ipp))+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]* sd.ipp.renew +fe4$coef[7]* sd.ipp.eGDP +fe4$coef[8]*mean(POLICYII$gasPrice))

mean.ipp;sd.ipp
100*(mean.ipp-sd.ipp)/mean.ipp

#effect of rps
mean.rps.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*mean(POLICYII$rpsprop)+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));mean.rps.renew
mean.rps<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*mean(POLICYII$rpsprop)+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]* mean.rps.renew +fe4$coef[7]*mean(POLICYII$ee) +fe4$coef[8]*mean(POLICYII$gasPrice))

sd.rps.renew<-pnorm(fixef(fe2)["CA"]+fixef(fe2)["CA"]*(mean(POLICYII$rpsprop)+sd(POLICYII$rpsprop))+fe2$coef[2]*0+fe2$coef[3]*mean(POLICYII$ipp)+fe2$coef[4]*0+fe2$coef[5]*mean(POLICYII$gdp)+fe2$coef[6]*mean(POLICYII$gasPrice));sd.rps.renew
sd.rps<-exp(fixef(fe4)["CA"]+fixef(fe4)["CA"]*(mean(POLICYII$rpsprop)+sd(POLICYII$rpsprop))+fe4$coef[2]*0+fe4$coef[3]*mean(POLICYII$ipp)+fe4$coef[4]*0+fe4$coef[5]*mean(POLICYII$gdp)+fe4$coef[6]*sd.rps.renew +fe4$coef[7]*mean(POLICYII$ee) +fe4$coef[8]*mean(POLICYII$gasPrice))

mean.rps;sd.rps
100*(mean.rps-sd.rps)/sd.rps
